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MRSURES  Oi  LACK  OF  FIT  FOR  RESPONSE  SURFACE  DESIGNS 
AND  PREDICTOR  VARIABLE  TRANSFORMATIONS+ 

G.  E.  P.  Box  and  N.  R.  Draper 


1 .  INTRODUCTION 


Consider  a  response  surface  study  in  which  a  polynomial  of  degree  d 

in  k  predictor  variables  . 5^.  is  used  to  represent  the  expected 

response  n  =  n( Ci » -  *  *  k ) •  Thus,  for  d  =  3,  we  have  a  third  order  model 


n  = 


k 

{  z 

i=l 


B!E->  + 
ri 


k 

{  I 


i=l  j>l 


k  k  k 

+  {  z  z  z 

i=l  j>i  Jl>i 


(1.1) 


Denote  the  observed  response  by  y  and  suppose  that  the  errors  eu  =  yu  “  Ou 

2 

are  independently  and  identically  normally  distributed  with  variance  a  .  If, 
as  is  usual,  we-  fit  models  of  this  kind  with  predictors  coded  in  "design  units" 


x,  .  <«,-«, „)'V 


(1.2) 


then  we  can  write  the  model  of  degree  d  =  3  in  the  form 


k  k  k 

H  =  3  +  {  Z  J3,x.)  +  {  Z  Z  B.  .x.x.) 
0  i=l  1  1  i=l  j>i  1J  1  3 

k  k  k 

+  {  Z  Z  Z  B-.^x.x.x^}. 

i=l  oii  &>i  J  J 


(1.3) 
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*odels  of  orders  d  =  1  and  d  =  2  are  obtained  by  omitting  the  appropriate 
Pracketed  terms  in  (1.1)  and  (1.3).  To  obtain  the  most  parsimonious  re¬ 
presentation,  (usually  corresponding  to  the  lowest  possible  value  of  d) 
it  may  be  necessary  to  transform  the  response  y  and/or  certain  of  the  pre¬ 
dictor  variables  £,. 

Now  whenever  we  fit  a  model  we  face  the  possibility  that  some  more 
complex  model  may  be  needed.  We  can  try  to  resolve  our  doubts  by  employing 
a  larger  design  which  makes  it  possible  to  fit  the  more  complex  models,  but 
then  similar  questions  arise  concerning  that  model  ,  and  so  on.  Clearly  we 
cannot  guard  against  all  possibilities.  A  practical  compromise  is 

1.  to  entertain  initially  at  each  stage  of  the  experimental  iteration, 
a  model  (containing,  say,  p  parameters)  which  it  is  hoped  will  be  adequate 

2.  to  employ  an  associated  design  which,  with  e  number  of  runs  only 
modestly  larger  than  p  provides  for  checks  sensitive  to  particularly 
feared  discrepancies. 

3.  if  such  discrepancies  occur,  to  consider  first  the  possibility  of 
their  elimination  by  transformation  (or  retransformation)  of  y  and/or 

^1  ’^2 '' ' '  ,f>k' 

4.  because  there  are  situations  where  a  more  complex  model  cannot  be 
avoided,  to  employ  design  arrangements  which  can  be  conveniently  augmented 
to  form  larger  designs  appropriate  for  fitting  and  checking  the  more  complex 
model . 
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Our  object  Is  to  consider  some  appropriate  checks  and  the  possible 
elimination  of  the  associated  lack  of  fit  by  power  transformations  of  the 
predictors  C»  for  first  and  second  order  polynomial  models  and  designs. 
Transformation  of  the  response  y  (Bartlett,  1947;  Box  and  Cox,  1964; 
Kruskal ,  1968;  Draper  and  Hunter,  1969)  will  not  be  considered  here. 


2.  FIRST  ORDER  MODELS  AND  DESIGNS 


2.1  Somp  first  order  designs  with  discrepancy  checks. 

Useful  first  order  response  surface  designs  are  the  two-level  factorials 
and  fractional  factorials  of  resolution  three  or  more  which  use  (respectively) 
all  or  some  of  the  2  runs  (+1  ,+l ,. . .  ,+l ) .  As  discussed,  for  example,  by 
Box  and  Wilson  (1951),  fractions  can  be  chosen  so  that  checks  are  associated 
with  the  residual  degrees  of  freedom  containing  feared  interactions.  Alter¬ 
natively,  or  in  addition,  (see  De  Baun,  1956),  by  adding  nQ  center  points  to 
any  such  design,  a  contrast  between  the  average  response  yCQ  at  the  center 
and  the  average  response  yc  at  the  factorial  points  is 

c2  =  *c  '  *co 

with  (2.1) 

k 

E(cJ  =  I  B... 

1  i=l  11 

Thus,  in  the  common  situation  where  the  B^  are  either  of  the  same  sign  or 
are  near  to  zero,  provides  an  overall  check  for  curvature  of  second  order. 

2 . 2  Can  we  use  a  first  order  model  in  transformed  predictor  variables? 

When  a  curvilinear  response  relationship  exists  whicn  is  monotonic  in 

the  predictor  variables  over  the  current  region  of  interest,  it  may  be  possible 
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^1  ^2 

to  use  a  first  order  model  in  which  power  transformations  Cj  >•••»£(< 
are  applied  to  the  C's. 

Assume  that,  at  worst,  the  response  may  be  represented  by  a  second 
degree  polynomial  in  transformed  variables,  namely  by 


k  X .  k 

n(C,X)  =  60  ♦  l  ^  ^ 


k 

E 


■  xi  xj 


(2.2) 


Then  a  first  degree  polynomial  model  will  be  appropriate  if  the  X^  may  be 

chosen  so  that  0- .  =  0  for  all  i  and  j.  In  Appendix  C,  we  show  that  this 
^  J 

requires  that 


where 


ni  j  =  0.  1  +  j 

(2.31 

6i(l-X.)ni  =  0,  i  =  l,2,...,k. 

(2.4) 

ni  = 


r  i 

r  .2  i 

an 

a  n 

_3x*. 

’  nij 

x.=0 

_3Vxj  J 

xr° 


(2.5) 


and  where 


(2.6) 
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Now  suppose  that  a  second  order  model  of  the  form  of  Eq.  (1.3)  with 
d  =  2  has  been  fitted  to  the  data  from  an  appropriate  design.  Then  we  could 
approximate  the  derivatives  of  Eq.  (2.5)  by 


ni  =  Vsi 


nii 


2bii/s;, 


(2.7) 


Thus  (i)  the  possibility  of  a  first  order  representation  in  transformed 
Ai 

variables  is  contra-indicated  (see  Eq.  (2.3))  if  one  or  more  interaction 
1  is  or 

estimates  b..,  i  4  are  significantly  different  from  zero,  and  (ii)  supposing 

1J  A 

such  a  transformation  to  be  possible,  the  appropriate  transformation  parameters 

★ 

are  roughly  estimated  by 


X,  =  1 


+  2bit  /<W* 


i  =  1 ,2,. . .  ,k. 


(2.8) 


More  precise  estimates  can  be  found  by  application  of  standard  nonlinear 

least  squares,  fitting  the  model  of  Eq.  (2.2)  with  8..  =  0  directly  to  the 

'  0 


data. 
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3.  SECOND  ORDER  MODELS  AND  DESIGNS 

3.1  Some  second  order  designs  with  discrepancy  checks. 

A  useful  class  of  second  order  designs  (see  Box  and  Wilson,  1951) 
appropriate  for  fitting  Eq.  (1.3)  with  d  =  2  consists  of  the  central  com¬ 
posite  arrangements  in  which  a  "cube",  consisting  of  a  two-level  factorial 
with  coded  points  (+1  ,+l , . . .  ,+l )  or  a  fraction  of  resolution  R  ^  5  is 
augmented  by  an  added  "star",  with  axial  points  at  coded  distance  a,  and 

by  nQ  added  center  points  at  (0,0 . 0).  More  generally,  both  the  cube  and 

the  star  might  be  replicated.  A  simple  example  of  a  design  of  this  type  for 
k  =  2  is  defined  by  the  columns  headed  x-j  and  x?  in  Table  1.  (We  use  nc 
and  n$  for  the  number  of  cube  and  star  points  and  nCQ  and  n$Q  for  the  number 

of  center  points  associated  with  them,  respectively.  Thus  n  =  n  +  n 

0  co  so  ' 

Also  shown  are  some  manufactured  data  whose  mode  of  generation  is  discussed 
in  Appendix  A. 

The  fitted  least  squares  second  degree  equation  is 

y  =  30.59  -  4.22x1  -  5.91x2  -  1.66x2  -  1.44x2  -  3. 41x^2  (3.1) 

+0.25  +0.17  +0.17  +0.14  +0.14  +0.24 

where  +  limits  beneath  each  estimated  coefficient  indicate  estimated 
standard  errors,  using  the  pure  error  estimate  s2  =  0.457  to  estimate  a2. 

An  associated  analysis  of  variance  table  is  shown  as  Table  2. 


Table  1.  A  composite  design  for  k  =  2  predictor  variables  and  its  associated 

estimator  columns;  n„  =  8,  n„„  s  1,  n  *  4,  n  =  5,  a  =  2. 

c  CO  s  so 
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Table  2.  Analysis  of  variance  associated  with  the  second  order  model  and 
its  checks, for  the  data  of  Table  1. 


Source 


df 


SS 


MS 


Mean 

Blocks 

First  order  extra 
Second  order  extra 
Lack  f  b. 


1 

1 

2 

3 

1 

1 

1 


16.74 
64.8'3<  173.17 
.14 


i 


Pure  error 


8 


3.657 


0.457 


-10- 


Before  accepting  the  utility  of  the  fitted  equation  we  would  need  to 
be  reassured  on  two  questions: 

1.  Is  there  evidence  from  the  data  of  serious  lack  of  fit?  If  not, 

A 

2.  Is  the  change  in  y,  over  the  experimental  region  explored  by  the 
design,  large  enough  compared  with  the  standard  error  of  y  to  indicate  that 
the  response  surface  is  adequately  estimated? 

The  analysis  of  variance  of  Table  2  sheds  light  on  both  these  questions. 
Its  use  to  throw  light  on  the  second  was  stuoied  by  Box  and  Wetz  (1973); 
see  also  Box,  Hunter  and  Hunter  (1978,  p.  524)  and  Draper  and  Smith  (1981, 
pp.  129-133). 

Clearly,  for  this  example,  it  is  the  marked  lack  of  fit  of  the  second 
order  model  that  immediately  concerns  us.  In  particular,  it  is  natural  to 
be  concerned  with  the  possible  effects  of  third  order  terms.  Associated  with 
the  design  of  Table  1  are  four  possible  third  order  columns  namely  those 
formed  by  creating  entries  of  the  form 

(x1,x1x2),  (x2,x2Xj).  (3..) 

These  form  two  sets  of  two  items,  as  indicated  by  the  parentheses. 

Now  suppose  these  third  order  columns  are  orthogonal i zed  with  respect 

to  the  lower  order  X-vectors.  This  may  be  accomplished  by  regressing  them 

against  the  first  six  columns  and  taking  residuals  to  yield  columns  x^^(from  x 
2 

x]22(from  x^x^),  and  so  on.  Then 
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m 


=  -3x . . . 
i  JJ 


t  J 


(3.3) 


and  the  residual  vectors  are  confounded  in  two  sets  of  two.  Furthermore, 
the  columns  and  ^22  are  orth°9°na^  t0  eac^  other.  These  vectors,  re¬ 
duced  by  a  convenient  factor  of  1.5  to  show  their  somewhat  remarkable  basic 
form,  are  given  in  Table  1. 

Consider  now  the  column  x^^  in  relation  to  Figure  1,  which  shows  the 
projection  of  the  points  of  the  composite  design  onto  the  x^  axis.  Denoting 
the  average  of  the  responses  at  x^  =  -a,  -1,  1 ,  a  by  y  Q,  y  ^ ,  y-j ,  and  ya, 
respectively,  we  see  that  a  contrast  c.^  associated  with  x^  is 


'31 


=  36 


1  ,ya“y-a 
r  2a 


yry- 


(3.4'; 


where,  for  our  example,  a  =  2.  The  expression  in  the  parentheses  is  an 
estimate  of  the  difference  in  slope  of  the  two  chords  joining  points  equi¬ 
distant  from  the  design  center.  For  a  quadratic  response  curve  this  difference 
is  zero.  Thus  c.^  is  a  natural  measure  of  overall  non-quadraticity  in  the  x^ 
direction.  A  corresponding  measure  in  the  x^  direction  is,  of  course,  given 

by  —  x222y<^' 

The  corresponding  sums  of  squares  for  these  contrasts,  given  in  Table  2, 
indicate  a  highly  significant  lack  of  fit.  Corresponding  plots  of  the  re¬ 
siduals  against  x-j  and  against  ^2  show  a  characteristic  pattern.  A  line 
joining  residuals  for  observations  at  x.  =  a  and  x.  -  -a  slopes  up,  while 
the  tendency  of  the  remaining  residuals  is  down  as  x.  is  increased.  We 
return  to  discuss  these  data  later. 
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Figure  1.  Projection  of  design  points  on  the  Xj  axis  for  the  composite  design 
In  two  factors  given  in  Table  1.  The  contrast  c^.  is  an  estimate  of 
the  difference  between  the  slopes  of  the  two  chords. 
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3.2  General  formulas. 

In  general,  a  composite  design  contains 

(a)  A  "cube",  consisting  of  a  2  factorial,  or  a  2  H  fractional  factorial, 

made  up  of  points  of  the  type  (+,+l . +1),  of  resolution  R  >_  5  (Box  and 

Hunter,  1961)  replicated  f (_>1 )  times.  There  are  thus  nc  =  f2k"p  such  points 
(where  p  may  be  zero). 

(b)  A  "star",  that  is,  2k  points  (+a,0,0 . 0),  (0,+a,0,. . .  ,0) , . . . , 

(0,0,0,. .. ,*x)  on  the  predictor  variable  axes,  replicated  r  times,  so  that 
there  are  n$  =  2kr  points  in  all. 

(c)  Center  points  (0,0, ...,0),  nQ  in  number,  of  which  nCQ  are  in  cube 

blocks  and  n„„  in  star  blocks, 
so 

It  is  shown  in  Appendix  B  that,  for  any  such  design,  k  sets  of 

2 

columns  can  be  isolated  with  the  ith  set  containing  the  k  columns  x.x., 

*  J 

j  =  l,2,...,k.  This  ith  set  is  associated  with  a  single  vector  x^.j  which 
is  orthogonal  to  the  (k+l)(k+2)/2  columns  required  for  fitting  the  second 
degree  equation  and  is  also  orthogonal  to  the  (k-1)  similarly  constructed 
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C3i 


5il1? 

FTTi^iii 


a2-l 


y  •  -y 

•7ai  -ai 
2a 


yli'y-li 
- T 


1 

] 


(3.5) 


with  standard  deviation 


Ct  J-  t 

3  a  -I 


1-+  1 


1/2 


nc  2ra2j 


(3-6) 


Also 


e(c31)  *  s(i1  *  (i-c2)-1  s.  e(ii. 


1JJ 


(3.7) 


and  the  contribution  to  the  lack  of  fit  sum  of  squares  is 


s;(c31)  =  (o2-!)2^/^*^) 


"3i  Ln 


c  2ra 


(3.6.) 


Note  that  even  if  E(c3i.)  =  0,  this  does  not  necessarily  mean  there  are 

no  cubic  coefficients.  A  combination  of  non-zero  8 . . i  and  g.  . .  could  occur 

2  -1 

for  which  g^..  +  (1-a  )  Ig..^  =  1S»  °f  course,  impossible  to  guard 

against  every  such  possibility  unless  the  full  cubic  model  is  fitted. 


3.3  Can  we  use  a  second  order  model  in  transformed  predictor  variables? 


In  some  instances,  lack  of  fit  of  a  second  order  model,  revealed  by 
significant  curvature  contrasts  of  the  kind  just  described,  might  be  removed 
by  transformations  of  some  or  all  of  the  predictors.  A  model  of  this  kind 
contains  many  fewer  parameters  than  a  full  third  order  model  and  is  much 
easier  to  analyze  and  interpret.  In  order  to  determine  conditions  that  must 
be  satisfied  to  make  it  possible  to  remove  lack  of  fit  in  this  way,  and  the 
part  that  the  curvature  contrasts  play  in  this,  we  suppose  that,  at  worst, 
the  response  function  may  be  represented  by  a  third  order  model  in  the  trans 
formed  predictor  variables.  In  Appendix  C  we  show  that  the  conditions 

S 

that  must  then  apply  if  all  third  order  coefficients  of  the  transformed 
are  to  be  zero  are 


nij£  =  °»  all  i  +  J  +  A-  1.2, —  ,k; 


(3.9) 


*  5f('-Xj,nij  ■  °-  =  . ki 


(3.10) 


niii  +  36i(l-xi)nii  +  sfO-x.Kl^A.)^.  =  0,  i  =  1,2,. ..,k. 


(3.11) 


An  important  conclusion  is  thus  the  following.  The  possibility  of  second 
order  representation  in  the  transformed  variables  is  contra-indicated  if  one 
(or  more)  interaction  estimates  b4  is  (are)  non-zero. 
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In  practice,  the  estimation  of  the  transformation  (when  not  contra¬ 
indicated)  is  best  done  using  nonlinear  least  squares  directly  on  the  model 
of  Eq.  (2.2);  however  some  interesting  light  on  how  the  curvature  measures 
c3i  relate  to  these  transformations  is  obtained  by  considering  how  Eqs.  (3.10) 
and  (3.11)  could  be  used  to  obtain  estimates  of  the  X^ . 

A  composite  design  does  not  permit  all  third  order  terms  in  Eq.  (1.3) 
to  be  separately  estimated.  Suppose,  however,  that  a  second  order  model 
augmented  with  only  cubic  terms 

eillxl  +  e222x2  +  "•  +  6kkkxk’ 


was  fitted.  If  the  response  n  could  be  represented  by  the  third  order  mcdel 
of  Eq.  (1.3),  the  estimates  and  b^  obtained  from  the  composite  design 
would  have  expectations 


E(bi)  =  ni  -  ^(l-a2)'1 

E(biii)  =  bui + 


ni  j  j 


'JJ 


(3.13) 

(3.14) 


If  now  b.  and  b^.  are  used  as  estimates  of  the  quantities  shown  as  their 
expectations  then,  after  appropriate  substitutions  have  been  made  in  Eqs. 
(3.9)-(3.11)  we  obtain  the  following  k  equations  for  the  X.. .  (In  these 


equations ,  bi ^  =  c3i . ) 
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U2, 


b-i.  +  6i(l-X.)bii  +^(1-Xi)(l-2Xi)b. 


l/i  ~2N -1 r,  1  2,2 


+  ^l-at)',{l-gat5^(l-2X.)(l-X.)  E  6.(  1  -X  . )b .  .}  =  0 
c  0  1  1  *  J  J  « J 


i  =  1 ,2,.. . ,k. 


(3.15) 


These  equations  can  be  solved  iteratively.  Guessed  values  for  the  Xi 
are  first  substituted  in  the  grouping  ( 1 -X^ ) ( 1 -2X . )  wherever  it  occurs  and  the 
resulting  linear  equations  solved  to  provide  improved  estimates  for  a  second 
iteration, and  sc  on. 

A 

For  the  example  data,  this  procedure  converges  to  the  values  X1  =  -0.23, 

^  A  /V 

X2  =  -0.93.  Those  may  be  compared  with  the  values  X1  =  0.09,  X2  =  -0.82 
provided  by  nonlinear  least  squares  (these  are  maximum  likelihood  estimates 
under  the  standard  normal  error  assumptions)  and  with  X^  =  0,  X2  =  -1,  the 
values  used  to  generate  the  data;  see  Appendix  A. 

An  analysis  of  variance  for  the  transformed  data  is  shown  in  Table  3 


where,  as  anticipated,  no  lack  of  fit  appears.  (See  also  Section  3.4  for  details.) 
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Table  3.  Analysis  of  variance  for  second  order  model  in  predictor 
variables  £n£^  and  • 


Source 

df 

SS 

MS 

Mean 

1 

13,938.934 

Blocks 

1 

7.347 

7.347 

First  order  extra 

2 

552.713 

276.357 

535.57 

Second  order  extra 

3 

565.238 

188.413 

365.14 

Lack  1 

\  CC 

f  1 

f  1.396 

f  1.396 

f  3.05 

< 

3  < 

2.021 

0.674 

1.47  j 

of  fit  i 

.Third  order 

L  2 

L  0.625 

L  0.313 

L  0.68 

Pure  error 

8 

3.657 

0.457 

Total 

18 

15,069.910 

3.4.  A  curvature  contrast 

Consider  the  overall  curvature  measure  C2  of  Eq.  (2.1)  used  to  check  the 
first  order  model.  When,  as  in  the  design  of  Table  1,  center  points  are  avail 
able  both  in  the  factorial  block(s)  and  in  the  star  block(s)  several  (two  for 
our  example)  such  measures  are  available.  Consider,  specifically,  the  two 
block  case  for  a  moment.  If  the  average  response  at  the  center  of  the  star 
is  y$0  and  the  average  over  all  the  star  points  is  y  ,  then  the  contrast 

c2  =  T  <V*so> 

a 
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has  expectation 


k 

E(cA)  =  Z  B . • 
6  1=1  11 


Thus  the  statistic 


cZ-Qi  V»c.  -  7 


which  is  the  difference  of  the  two  measures  of  overall  curvature,  should  be 
zero  if  the  assumptions  made  about  the  model  being  quadratic  are  true. 

From  Figure  1,  we  see  that  the  curvature  measure  c^  associated  with  tne 
cube  (open  circles)  is  contrasted  with  c£  associated  with  the  star  (block  dots). 

In  general  the  distance  from  the  center  of  the  design  to  the  cube  points  is 

1/2  1/2 
k  and  that  for  the  star  points  is  a.  When,  as  in  our  example,  k  and  a 

are  different,  a  significant  value  of  ^2~c2  cou^ci  indicate  (for  example)  a 

symmetric  departure  from  quadratic  fall -off  on  each  side  of  the  maximum,  such 

as  we  see  (for  example)  in  a  normal  distribution  curve. 

In  general,  for  two  blocks,  the  standard  deviation  for  c^-d,  is  given  by 


o  ,  .  {J- ♦  J- ♦  4  [  '  ♦ 

c2 -°i  "c  "co  c4  l2kr  "so 


and  the  associated  sum  of  squares  for  the  analysis  of  variance  table  entry 
of  Table  2  is  obtained  from 
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SS^-c^)  =  /{  +  +  4  +  4  } 

c  nc  n  co  2rcT  a  r> 

so 


For  our  example  we  find 


c2-c£  =  1.3925  +  0.7520 


with  associated  sum  of  squares  1.567  as  shown  in  Table  2.  There  is  clearly 
no  evidence  of  this  sort  of  lack  of  fit.  When  transformed  predictors  are  used, 
again  no  lack  of  fit  of  this  kind  is  evident,  as  we  see  from  Table  3. 


3.5.  Interaction  with  blocks? 

When  composite  designs  are  run  in  blocks,  and  if  we  allow  the  possibility 
that  effects  from  the  predictor  variables  could  interact  with  blocks,  then 
the  various  measures  of  lack  of  fit  would  be  confounded  with  block-effect 
interactions.  Although  such  contingencies  must  always  be  borne  in  mind,  it 
should  be  remembered  that  these  particular  block-effect  interactions  are  no 
more  likely  than  any  others. 


L 
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SUMMARY 


In  the  traditional  analysis  of  second  order  response  surface  designs,  a 
number  of  degrees  of  freedom  are  usually  consigned  to  "lack  of  fit"  and  the 
corresponding  sum  of  squares  is  used  to  test  for  overall  lack  of  fit.  Some 
split-up  of  lack  of  fit  has  been  previously  discussed  (Draper  and  Herzberg, 
1971),  but  a  much  more  ambitious  and  detailed  division  is  described  here.  We 
show  that  it  is  possible  to  check  for  cubic  lack  of  fit  in  the  k  axial 
directions  and  if  it  exists,  not  only  to  chect  if  it  is  possible  to  eliminate 
it  by  a  power  transformation  in  the  predictors,  but  also  actually  to  estimate 
the  powers  needed  to  effect  the  transformation.  The  theory  is  derived  in 
Appendix  C,  and  a  worked  two-factor  example  shows  how  to  carry  out  the  cal¬ 
culations. 

It  is  also  shown  how  a  certain  curvature  contrast  can  be  used  to  check 
overall  quadratic  fall-off  away  from  the  maximum  of  a  response  surface. 

Simpler  but  similar  considerations,  one  degree  down,  apply  to  the  first 
order  model,  and  appropriate  formulas  for  these  are  derived  as  a  special  case 
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Appendix  A.  Generation  of  example  data. 

Figure  AT  is  taken  from  the  manuscript  of  a  book  on  response  surfaces 
by  G.E.P.  Box,  N.R.  Draper,  and  J.S.  Hunter,  in  preparation.  Figure  Al(b) 
shows  a  quadratic  response  function  with  a  simple  maximum  in  variables 
Un^,  lOOC^)  •  This  figure  is  redrawn  in  the  metrics  in  Figure 

Al(a).  In  the  representation,  the  doubled  cube  plus  star  plus  center 

points  design  of  Table  1  is  indicated  by  the  positions  of  the  dots.  Response 
values  were  calculated  at  these  points,  and  random  error  added  to  give  the 
y-values  in  Table  1.  For  these  generated  data,  Sjq  =  2.5,  E^g  =  12.5, 

S-j  =  0.75,  and  S£  -  3.75. 
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Appendix  B.  The  third  order  columns  of  the  X-matrix  for  a  composite  design. 

We  here  prove  the  results  summarized  in  Section  3.2.  The  full  cubic 
model  in  k  variables  Xj^.-.-.x^  is  given  by  Eq.  (3.2).  The  form  of  the 
X  matrix  in  the  regression  model  y  =  XB  +  e  when  the  design  consists  of  f 
"cubes"  plus  r  "stars"  plus  nQ  center  points  is  as  shown  in  Table  B1 . 

We  can  denote  columns  by  placing  square  brackets  around  the  column  head;  for 
example  [x^]  will  denote  the  x-j  column,  and  so  on.  We  write  nc  =  f2  p  for 
the  number  of  cube  points. 

All  of  the  cubic  columns  are  orthogonal  to  all  of  the  other  columns 

3  2 

with  the  following  exceptions;  [x?]  is  not  orthogonal  to  [x^],  nor  to  [x^x^]; 

2  3  2 

[x.x^]  is  not  orthogonal  to  [x.],  nor  to  [x.],  nor  to  [x.x  ].  The  first  step 

1  J  1  *  1 

3  2 

is  to  regress  the  [x^  and  [x..Xj]  vectors  on  the  [x^  and  take  residuals. 
Because  the  columns  involved  are  orthogonal  to  [xQ],  no  adjustment  for  means 

O 

is  needed.  We  denote  the  "cube  portion"  of  the  [x.]  and  [x^xj]  vectors  by  c^ , 
as  indicated  in  the  table.  These  two  sets  of  residuals  are,  where  the  prime 
denotes  transpose, 

[x11i]  =  Cxf]  -  {[xi]1[xj]/[xi]*[xi]}[xi]  (Al) 

and 

[Xijj]  =  O.jX'j]  -  {[xi]'[xix2]/[xi]I[xi]}[xi]  (A2) 

p 

both  of  which  reduce  to  multiples  (1-a  )m  and  m,  respectively,  where 
2  2 

m  =  2ra  /(nc+2ra  ),  of  the  same  vector.  For  example, 


Table  Bl.  The  X  Matrix  Tor  n  Cubic  Model  in  k  Predictors 
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[x,,,]'  =  [cj,d,-d,d,-d,...,d,-d,0,0,...,0](l-a2 )m, 


(A3) 


where 


d  =  nc/(2ra),  (A4) 

and  where  there  are  r  sets  of  ( d * -d ' s )  in  the  vector.  In  general,  for  [x.^]' 
will  be  replaced  by  c^‘  and  the  position  of  the  +d's  will  correspond  to  those 
of  the  +a's  in  the  corresponding  [x^]'  vector.  Note  that,  because 
c.'c,  =  0,  i  f  j,  it  is  obvious  that  Cx .  -  -  ]  and  [x...]  are  orthogonal. 

1  J  III  J  J  J 

It  follows  that  the  k  cubic  coefficients  B^,  B^j  ( j + i  ,j=l  ,2,. . .  ,k, 
otherwise)  cannot  be  estimated  individually  but  only  in  linear  combination, 
and  that  an  appropriate  normalized  estimating  constrast  for  this  is 


{$i^i  +  d(-ryai+ry.ai)>/(nc(1-a  )> 


(A5) 


where  y^  is  the  portion  of  y  corresponding  to  the  cube  part  of  the  design, 
and  y^,  y_a^  are,  respectively  the  averages  of  observations  taken  at  the  a 
and  -a  axial  points  on  the  x..  axis.  If  we  similarly  denote  by  y^  and  y_li. 
the  averages  of  the  nc/2  observations  in  y^  corresponding  to  1  and  -1  in  c^ , 
respectively,  it  follows  quickly  that  =  c^.  where  c^  is  given  in  Eq.  (3.5) 


The  expected  value  is 
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E(c3t>  ' 


where  X  is  as  in  Table  B1  and  the  coefficients  of  B  correspond  to  the 

columns  in  the  obvious  manner.  Because  [x^..]  is  orthogonal  to  all  columns 

3  2 

of  X  except  the  [x.]  and  [x.x.]  columns,  Eq.  (3.7)  emerges  almost  immediately. 

^  '  *  J 


Appendix  C.  Conditions  for  efficacy  of  transformations  on  the  predictor 


variables. 

Suppose  that  an  experimental  design  is  run  in  the  k  coded  variables 


X,  -  (5r<io>'V 


1  =  12  k 

I  I  )L  9  •  •  •  9«S  9 


(Cl) 


in  a  situation  where  the  underlying  response  function  can  be  approximated 

Xi 

by  a  second  degree  polynomial  F{(£..  )}  in  the  transformed  original  variables 

i  - 

Ai 

Ci  .  Thus 

k  A.  k  k  A.  A. 

F(5,X)  -  6  +  I  M*  +  E  E  (C2) 

c  i=l  1  1  i=l  j>1  1J  1  3 


Assume  all  A^  j  0.  (The  case  when  any  Ai  *  0  can  be  handled  as  a  limiting 
case  using  the  fact  (see,  for  example,  Box  and  Cox,  1964)  that  (£X-1)/A 
tends  to  X.nA  as  A  tends  to  zero.) 

If  (C2)  is  to  be  a  suitable  representation,  then  all  third  derivatives 
A. 

with  respect  to  the  C..1  must  vanish  identically.  Note  that 


3F  9F  3Ci 
Ai  -  3^  A  . 
3C- 


1  -A . 

=  f,U<  Va.) 


(C3) 


say,  where  F..  =  3F/3C^ *  Moreover,  because  of  Eq.  (Cl), 

i 
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_  3F  3F  9Ci  _  -  c 
ni  "  3x.  "  3^i  3x1  "  ii 


(C4) 


so  that 


3F 


1-X, 


35, 


(C5) 


The  obvious  extensions  of  these  results  also  follow  for  the  higher  derivatives 
which  involve  terms  of  the  form 


32F 


*i  j  Sx^Xj 


and  n. 


33F 


ijJL  Sx^SXjSx^ 


(C6) 


If  we  now  carry  out  the  appropriate  differentiations,  we  obtain 

l-*i 

3F  _5i 


35 


x.  -  S.X.  ni 


(C7) 


32F 


1-X  1-X 

£.  V.  J 

'  JL 


^ '  Wj 


(C8) 


d  r 


(V) 


( C9) 


92F 


A .  2  ^ 

9  (CjJ) 


1-A.  1-2A 

s.j 

VjVj 


=  0 


(CIO) 


33F 


V> 


1-3A. 

*i 

5 

Vi 


(57)^111  +  I7-  (1-^i)nii  -  (l-Ai)(l-2Xi)ni 


=  0 
(Cll) 


33F 


A  A  A 


1-A.  1-A,  1-A 

Ci  *1  \ 

ViViYs, 

1  J  ~  *  J  ** 


l 


'ij£ 


=  0. 


(C12) 


Eqs.  (CIO)  -  (C12)  are  exact  conditions  on  the  A's.  The  and 

also  involve  the  A's,  and  cannot  be  specifically  evaluated  if  the  A's  are  un¬ 
known.  If  estimates  of  the  n-,  and  n.-in  are  substituted,  however,  Eqs. 

-  1  1J  1JX, 

(CIO)  -  (C12)  can  be  solved  to  provide  estimates  of  the  A's. 

We  now  assume  that  the  response  surface  can  be  approximately  represented 


by  a  cubic  polynomial  in  the  coded  predictor  variables,  namely,  by 
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n  =  bo  +  bl*l  *  •••  *  Vk 

+  bllxl  +  •••  +  bkkxk 
+  b12x]x2  +  ...  +  bk_1>|cxklxk 

+  b^x^  +  b]22x1x2  +  ...  (C13) 

3  2 

+  b222x2  +  b-j  1 2xi  x2  +  •  •  • 

+  bkkkxk  +  bllkxlxk  +  "* 

+  b123x1x2x3  +  ... 


We  can  estimate  the  n- .  n.-.-,  and  n-  •  »  by  the  corresponding  derivatives  of 

1  I J  1 J  ** 

A 

n  evaluated  at  the  center  of  the  design,  that  is,  at  xi  =  0.  In  general  then, 


ni  =  bi  ’  nii=2bii’  niii=6biii’ 


(C14; 


nij  =  bi.i’ 


nijj  =  2bijj’  niji  “  bijJt* 


and  we  would  substitute  these  in  (CIO)  -  (C12),  at  the  same  time  setting 

£.  =  C-_»  its  value  when  x.  =  0.  Note  that,  from  Eq.  (Cl 2),  we  require  that 

1  10  1  A. 

n. . .  =  0,  which  implies  that  the  assumed  transformations  C-1  are  not  suitable 
1 * 

representations  unless  all  three-factor  interactions  are  zero.  In  practice, 
then,  we  would  want  b.-0  to  be  small  and  nonsignificant  when  i  f  j  f  l  or 
else  we  have  a  clear  indication  that  the  will  not  provide  a  satisfactory 
second  order  representation  in  (C2).  If  this  aspect  is  satisfied,  however, 
we  now  solve  the  ^k(k+l)  +  k  =  ^k(k+3)  simultaneous  equations: 
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"ijj  *  {j(1-XJ>n1j  "  °-  . k* 


niii  +  36i(l-Ai)nii  +  at(l-Ai)(l-2xi)n1  =  0, 


i  =  1,2 . k  (C15) 


where  6.  =  and  we  have  divided  through  Eqs.  (CIO) 
factors  assumed  to  be  non-zero,  namely  ,  using 
Eq.  (C14) . 


and  (Cl  1 )  by 

A 

the  q  values  from 


Composite  Designs 

An  additional  complication  arises  with  composite  designs.  For  such  de¬ 
signs,  we  cannot  estimate  all  the  third  order  coefficients  individually,  a:; 
we  have  explained  in  Appendix  B.  This  means  that,  while  Eqs.  (Cl 5)  are  still 

valid,  the  values  in  (C14)  cannot  be  used.  For  a  composite  design,  the  column 
3  2 

vectors  x..,  x. ,  and  x.^  are  linearly  dependent  via 


x.x2.  =  (-a^.+x^/O-a2). 


(C16) 


Thus  in  the  general  cubic  model,  we  cannot  estimate  all  k+k+k ( k- 1 )  =  k(k+l) 
coefficients  in  the  terms 


6i*i 


8. • -x? 
m  i 


k 

l 

J'fi 


B.  •  .x.x. 
UJ  i  J 


(C17) 


but  only  the  2k  coefficients  of 


2  k  ,  k  3 

18.-  — Z  +  {SiH  +  “S'  2  6iii}xi 

1  1 -a  jfi  TJJ  1  111  1 -or  jfi  1JJ 


1-^  j+1 


(CIS) 


It  follows  that,  in  the  model  so  reduced. 


a2  1  k 

b.  estimates  n, - 5-3-  2  ^ 

1  1  1-a2  2  jfi  1JJ 


(Cl  9) 


while 


bu(  estimates  ^  r,... 


(C20i 


(By  examining  Eqs.  (3.7),  (C14) ,  and  the  second  portion  of  Eq.  (C18),  we 

infer  that  b.^.  =  c^ . )  Alternatively,  if  the  model  is  fitted  using  and 

3  3 

the  "orthogenalized  x^",  namely  x^.  =  x^  -  4>x. ,  the  terms  of  the  model  are 

<ei  +  *eiii  *ej,  6ioj,xi +  ,siii  *  A  j,  Bijj,ui  -  »xi>  <c2" 


where 


8  =  nc^nc+2m2'  ’  ^  =  (nc+Zm4)/(nc+2ra2)  •  (C22) 


In  this  form,  we  have  that 
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b*  estimates 


6^niii 


1  ^ 

1  0  Z 

2 


Ijj 


(C23) 


and 


b*..  estimates  ^n...  + 

111  Dili 


l 


l  k 
l-a  2  jf i 


nijj 


(C24) 


3 

where  b|  is  the  estimated  coefficient  of  x..  and  b*..  is  that  for  (x..  -  yx.). 
(Note  that  bt..  =  b...  =  c~. ,  but  that  b*  =  b.  +  *Pb . . . . )  We  now  describe 

ill  I  |  1  O  I  i  I  III 

how  this  affects  the  estimation  of  the  A...  From  Eqs.  (CIO)  and  (Cll),  and 
setting  £.  =  £.  we  have 


nijj  *  6j(,-xj)nij  =  0 


(C25) 


and 

ni11  +  36i(l-Xi)nii  +  ^(l-A-W^A.)^  =  o  (C26) 

2 

We  now  combine  ^  times  (C26)  with  { - — ~ - - — x—  6.(1  -2X - )(1-X.)>  times 

6  2(l-a  )  12(1 -a2)  1  1  1 

(C25)  summed  over  j  f  i  to  give,  for  i  =  l,2,...,k, 


r 
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i 

i 


+  *  ii2  jti 

-  A?  ].  n^.j)  (C27) 

1  -Ot  jfl 

z  6iO-xi)  n.-i  =  0. 

1  l-a  b  1  1  1  jfi  J  J  1J 

O 

Thus,  if  the  third  order  model  is  fitted  in  terms  of  and  x":  (rather  than  x. 
3 

and  (xi  -  ^x. )  as  described  above)  we  can  substitute  appropriate  estimates 
to  give  the  k  simultaneous  equations  for  i  =  l,2,...,k. 


1,2/ 


biii  ♦V^Vbii  ^o-VO^Vbi 


1  1 


12 ,2 


+  ^-Vn-^x^Jd^XJd-X.)}  Z  6  .  ( 1  -X  . )  b . . 
1  iV  51  1  1  j+i  J  J  1J 


=  0 


(C28) 


which  can  be  solved  for  X-j  .Xg,. . .  ,X^.  These  awkward  equations  can  be  difficult 
to  solve  unless  some  care  is  applied.  We  suggest  an  iterative  procedure  in 
which  rough  estimates  of  X^  are  used  in  the  grouping  Qi  =  (1-X.)(1-2X^) 
which  occurs  in  two  positions  in  Eq.  (C28).  The  resulting  linear  equations 
in  6i  =  1  -X..  are  straightforward  to  solve,  and  the  results  are  used  in  the 
grouping  (1-X^)(1-2X. )  for  a  second  iteration  and  so  on,  until  convergence  is 


-36- 


achieved.  To  aid  convergence,  each  new  iteration  can  be  started  from  the 
midpoint  of  the  old  and  new  values,  if  desired. 


Alternative  Cubic  Case 


3 

If  the  cubic  is  fitted  with  estimated  terms  b*x.  +  bt„ (x^x..)  as 
described  above,  we  obtain  b^  and  b^ from 


b. 


bt  -  ipb* . .  and 
i  m 


b...  =  b*.., 
in  m 


(C29) 


Related  Work 


The  expression  in  (C20)  can  be  alternatively  written  in  the  form 


0. . .  +  (l-a2)_1ZB _ 

■m  v  '  i.ij 


Draper  and  Herzberg  (1971,  p.  236)  show  that  the  sum  of  squares  of  these 
quantities  for  i  =  l,2,...,k  occurs  in  the  expected  value  of  a  general  measure 
of  lack  of  fit  Lj.  Thus,  the  c^  contrasts  essentially  provide  a  split-up  of 
which  permits  a  more  detailed  and  sensitive  analysis.  The  remaining  degrees 
of  freedom  pertaining  to  L-j  can  be  attributed  to  other  contrasts  as  already 


described. 
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The  First  Oegree  Case 


We  return  to  the  beginning  of  this  Appendix  to  examine  the  simpler  case 

when  the  underlying  response  surface  can  be  approximated  by  a  first  degree 

X. 

polynomial  in  the  transformed  variables  ^ 1 .  In  such  a  case 

1.  All  8'  =  0,  in  (C2). 

X. 

2.  All  second  derivatives  with  respect  to  must  vanish  identically. 

This  condition  provides,  from  Eqs.  (C8)  and  (09),  the  equations  (2.3)  and  (2.4). 
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